This invention relates generally to methods, apparatus and software for image reconstruction.
Traditionally, images have been reconstructed from computed tomography (CT) data using so-called direct reconstruction algorithms such as filtered back projection (FBP) or convolution back projection (CBP). Recently, iterative reconstruction (IR) algorithms have been introduced for the reconstruction of CT images. One advantage of IR algorithms is that they can more accurately model the measurements obtained from real CT systems. This is particularly true for helical CT systems with multi-slice detectors because these systems produce projection measurements that pass obliquely through the 2-D reconstructed image planes. By more accurately modeling these projections, IR algorithms can produce reconstructions with higher quality, lower noise, and fewer artifacts.
For example, consider a helical scan CT system. The 3-D volume to be reconstructed can be represent by an array of N discrete voxels xi where i is the index of the voxel 3-D position. The value xi may specify the unknown density of the voxel. Furthermore, let x=[x1,x2, . . . ,xN] be a vector containing the unknown density of each voxel in the reconstruction. So in this case, x represents the full 3-D reconstruction volume. During the CT scanning process, projections are measured for M different projections through the object. The different projections are typically measured for a wide variety of positions and angles through the object. The value of the integral jth projection through the object is denoted by ym and the vector of all measurements is denoted by y=[y1,y2, . . . yM].
The measurement ym of the CT scanning process can be modeled as being a noisy version of the forward projection of the 3-D volume x. So the mth projection can be analytically computed as Fm(x) where Fm(●) is the forward projection function. The difference between ym and Fm(x) is caused by noise and distortions when x is chosen to exactly match the true voxel densities of the object being imaged.
The mismatch between the true projection measurement ym and the analytical forward projection Fm(x) can be evaluated using a distance function of the form Dm(ym,Fm(x)) where the distance function Dm(●,●) penalizes mismatch between ym and Fm(x). The distance function Dm(●,●) should be large when ym and Fm(x) are very different and small when ym and Fm(x) are very similar. The distance function is not required to obey any specific metric properties such as the triangle inequality. In some cases, it can be chosen to have the form Dm(ym,Fm(x))=wm|ym−Fm(x)|2 where wm is a scalar weight. In some cases, it may be chosen to have the form Dm(ym,Fm(x))=p(ym|Fm(x)) where p(●|●) is a conditional probability function for a Poisson, Gaussian, or other distribution.
The objective of IR algorithms is to determine the unknown value of x by searching for the value of the vector x that best matches the measured data. Typically, this is done by minimizing a cost function of the form
                              x          ^                =                  arg          ⁢                                          ⁢                                    min              x                        ⁢                          {                                                ∑                                      m                    =                    0                                    M                                ⁢                                                                  ⁢                                                      D                    m                                    ⁡                                      (                                                                  y                        m                                            ,                                                                        F                          m                                                ⁡                                                  (                          x                          )                                                                                      )                                                              }                                                          (        1        )            where {circumflex over (x)} is the value of the variable x which achieves the minimum of the function. This cost function can be minimized in a variety of manners using optimization methods such as iterative coordinate descent, expectation maximization, conjugate gradient, or any number of alternative techniques. In practice, the solution to (1) is often too noisy for clinical use.